'''
Created on Jun 15, 2010

@author: oabalbin
'''

import sys
import copy
import pickle
import numpy as np
from optparse import OptionParser
from datetime import datetime
from collections import defaultdict, deque

import amarillo.copa.copa_analysis as ca
import amarillo.copa.subset_genes as sg
import amarillo.parser.parse_gene_lists as pg
import RNAseq.array as seqarray

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    optionparser.add_option("-s", "--normalsFile", dest="normalsFile",
                            help="annotation file with the normal samples")
    optionparser.add_option("-t", "--tissueName", dest="tissueName",
                            help="name of the tissue analyzed")


 

    
    (options, args) = optionparser.parse_args()
    
    
    t = datetime.now()
    tstamp = t.strftime("%Y_%m_%d_%H_%M")
    
    outfiles = ['outfile_outliers','outfile_outliers_names','outfile_outliers_list','outfile_outliers_matrix','outfile_copa_matrix',
                'outfile_foldchange_matrix','outfile_foldchangesinsample']
    time_stamp_outfiles=[]
    for name in outfiles: 
        time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.tissueName+'_'+name+'.txt')
        
    thisfile = open(options.annotFile)
    outfile_outliers= open(time_stamp_outfiles[0],'w')
    outfile_outliersinsample= open(time_stamp_outfiles[1],'w')
    #outfile_outliers_list = open(time_stamp_outfiles[2],'w')
    outfile_outliers_matrix = open(time_stamp_outfiles[3],'w')
    array_dump_file = options.annotFile + '_dump_file'
    outfile_copa_matrix = open(time_stamp_outfiles[4],'w')
    outfile_foldchange_matrix = open(time_stamp_outfiles[5],'w')
    outfile_foldchangesinsample = open(time_stamp_outfiles[6],'w')
    
    #gene_list_folder = '/data/metadata/panther_gene_list/'
    gene_list_folder ='/data/metadata/my_kinases_list/'
    #gene_list_folder = '/data/metadata/my_gene_lists/'
    #samples_list_file='/data/projects/outliers/prostate_tissue_samples.dat'
    samples_list_file='/data/projects/prelims/lung_cell_lines.txt'
    #samples_list_file='/data/projects/outliers/lung_cell_lines.txt'
    #samples_list_file='/data/projects/outliers/breast_celllines_chandan2.dat'
    #gene_list_file = '/data/projects/outliers/kinases_list_science2002.dat'
    ####################
    ##### Parameters
    
    normal_samples_list = ca.list_of_names(open(options.normalsFile))
    
    #percentils_list = [50]
    percentils_list = [75,80,85,90, 95] #[80,85,90,95,98] #[75,80,85,90, 95] #[80,85,90,95,98] 
    method='mean' # median, mean, sum
    #outliers_cutoff=100 #100 #200
    nout2report=10
    expression_cutoff=0.0       # minimum median RPKM of expression
    base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                                # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
    median_shift=1.0
    sdv_factor=0.25             #0.25
    not_report_isoforms=True   # Report gene outlier isoforms
    
    ########### Program
           
    # create the meme array
    filename=options.annotFile
    matrixfilename=options.annotFile+'memmap'
    header_cols=2
    matobj = seqarray.read_matrix(filename, matrixfilename, header_cols)
    all_samples, all_genes, all_gene_intervals = \
    seqarray.create_matobj_dicts(matobj.cols,matobj.sample_names)
    
    #####
    
    # Get all the gene files in the folder
    thosefiles = sg.read_files_folder(gene_list_folder,'.txt')
    
    
    thissampleslist = ca.list_of_names(open(samples_list_file))
    sample_subset, sample_cols, sample_subset_rev = sg.get_sample_subset(thissampleslist, all_samples)
    #normal_sample_subset, normal_sample_cols, normal_sample_subset_rev = sg.get_sample_subset(normal_samples_list, all_samples)
    normal_sample_subset, normal_sample_cols, normal_sample_subset_rev = sg.get_sample_subset(normal_samples_list, sample_subset_rev)

    #cohort_outlier_profile=defaultdict(deque)
    all_outliers_intersection={}
    cohort_outlier_profile_dict={}
    all_gene_list_names=deque()
    all_unique_genes = []

    for gene_list_file in thosefiles:
        
        gene_list_name = gene_list_file.split('/')[-1]
        all_gene_list_names.append(gene_list_name)
        
        outfile_outliers_list = open(time_stamp_outfiles[2]+'_'+gene_list_name,'w')
        
        # subset of genes and samples to consider
        thisgenelist = ca.list_of_names(open(gene_list_file))
        #thisgenelist = pg.panther_gene_list(open(gene_list_file))
        # Only unique genes not considered previously in other lists
        thisgenelist = list(set(thisgenelist).difference(set(all_unique_genes)))
        [all_unique_genes.append(g) for g in thisgenelist]
        
        #thissampleslist = ca.list_of_names(open(samples_list_file))
        gene_subset, gene_rows, gene_subset_intervals = sg.get_gene_subset(thisgenelist, all_genes, all_gene_intervals)
        #sample_subset, sample_cols = sg.get_sample_subset(thissampleslist, all_samples)
        #normal_sample_subset, normal_sample_cols = sg.get_sample_subset(normal_samples_list, all_samples)
        
        # create a sub array  using the genelist and the samplelist
        medthr=1.5
        expression_matrix = matobj.exprs[gene_rows, :]
        expression_matrix = expression_matrix[:,sample_cols]
        '''
        indicator = np.array(range(len(expression_matrix)))
        print len(indicator)
        matrix_med = np.median(expression_matrix,axis=1)
        indicator = indicator[matrix_med >= medthr]
        print len(indicator)
        '''
        
        thisGeneArray = sg.create_thisGenearray(expression_matrix, gene_subset, sample_subset, \
                                                list(normal_sample_cols), gene_subset_intervals)
                
        thisGeneArray.copa_norm_across_samples(median_shift)
        thisGeneArray.copa_norm_in_sample()
        thisGeneArray.fold_change_accross_samples()
        #outliers_cutoff = thisGeneArray.expVal.shape[0]
        outliers_cutoff = np.ceil(0.01*len(gene_subset))
        outliers_cutoff = 50
        #print outliers_cutoff
        
        ############# Print expression matrix after Copa normalization
        ca.print_copa_matrix(thisGeneArray, outfile_copa_matrix)
        #ca.print_foldchange_matrix(thisGeneArray, outfile_foldchange_matrix)
        #ca.print_foldchange_outliers_by_samples(thisGeneArray,outfile_foldchangesinsample, nout2report)
        #sys.exit(0)
        ############ Run the outlier analysis for different percentiles
        
        gene_outliersample_dict, copa_rank_dict, copa_score_dict = ca.copa_at_different_percentiles(thisGeneArray,percentils_list)
        
        ############ Summarize the results from the different percentiles
        #_ETS
        final_list_genes, gene_copascore_values = ca.summarize_results_copa_percentiles_ETS(thisGeneArray,
                                                                                      gene_outliersample_dict, copa_rank_dict, 
                                                                                      copa_score_dict, percentils_list, sdv_factor, method)
        
        ############# Identify all outliers for each sample
        # sorted outliers to report and to matrices
        list_outliers_report, outliers_in_samples, outliers_in_samples_copascore = ca.identify_outliers_in_sample(thisGeneArray,final_list_genes, gene_outliersample_dict, 
                                                                                                                  gene_copascore_values, outliers_cutoff, not_report_isoforms)
        
        
        ### Print outlier list and outlier matrix. 
        ca.print_list_of_outliers(list_outliers_report, outfile_outliers_list)
        ca.print_matrix_of_outliers(thisGeneArray.get_list_of_samples(), list_outliers_report, 
                                    outliers_in_samples_copascore, outfile_outliers_matrix)
                
        
        ############# Print sample : outliers gene in that sample
        # Returns a dictionary with sample and the outlier in that sample for each sample.
        cohort_outlier_profile, outliers_sample_intersection =  ca.get_outliers_matrix_by_samples_mod(thisGeneArray, gene_copascore_values, 
                                                                                    outliers_in_samples, outfile_outliersinsample, nout2report)
        cohort_outlier_profile_dict[gene_list_name] = cohort_outlier_profile
        all_outliers_intersection[gene_list_name] = outliers_sample_intersection
        
        thisGeneArray =[]




outliers_intersection_dict=defaultdict(set)
#outfile_outliersinsample.write('sample'+'\t'+",".join(outliers_header).replace(',','\t')+'\n')
list_of_genes = []
list_of_genelistnames=[]
list_of_values= defaultdict(list)
for thisglname in all_gene_list_names:
    this_cohort_outlier_profile = cohort_outlier_profile_dict[thisglname]
    this_intersect_profile = all_outliers_intersection[thisglname]    
    
    tempgenes, tempnames=['z' for i in range(len(this_intersect_profile))],[]
    for g, i in this_intersect_profile.iteritems():
        tempgenes[i]=g
        print g,i     
    #print tempgenes
    list_of_genes.append(list(tempgenes))
    [tempnames.append(thisglname) for i in range(len(this_intersect_profile))]
    list_of_genelistnames.append(tempnames)
    
    
    for sample, outlier_profile in this_cohort_outlier_profile.iteritems():
        # Make a new array for registering outlier values
        thisarray =  np.zeros(len(this_intersect_profile))
        
        for pair in outlier_profile: 
                       
            name = pair[0]
            if name!='NAN':
                gen_ind = this_intersect_profile[name]
                thisarray[gen_ind]=pair[1]
            
        #print sample, map(str,list(thisarray))
        
        list_of_values[sample].append(list(thisarray))
        #outfile_outliersinsample.write(sample+'\t'+",".join(map(str,list(thisarray))).replace(',','\t')+'\n')    
        #outfile_outliersinsample.write('\n')
        
outfile_outliersinsample.write('sample'+'\t'+",".join(sum(list_of_genelistnames,[])).replace(',','\t')+'\n')
outfile_outliersinsample.write('sample'+'\t'+",".join(sum(list_of_genes,[])).replace(',','\t')+'\n')
for sample, profile_values in list_of_values.iteritems():
    #print sample, sum(profile_values,[])
    outfile_outliersinsample.write(sample+'\t'+",".join(map(str,sum(profile_values,[]))).replace(',','\t')+'\n')
    
    
    
    
    
    